home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / brent_cc.lha / brent_cc / zeroin.cc < prev   
C/C++ Source or Header  |  1993-08-08  |  5KB  |  143 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *              Numerical Math Package
  6.  *
  7.  *                Brent's root finder
  8.  *           obtains a zero of a function of one variable
  9.  *
  10.  * Synopsis
  11.  *    double zeroin(ax,bx,f,tol=EPSILON)
  12.  *    const double ax         The root is to be sought within
  13.  *    const double bx          the interval [ax,bx]
  14.  *    double (*f)(const double x)    Ptr to the function under 
  15.  *                    consideration    
  16.  *    const double tol        Acceptable tolerance for the root
  17.  *                    position. It is an optional parameter
  18.  *                    with default value EPSILON
  19.  *
  20.  *    Zeroin returns an approximate location for the root with accuracy
  21.  *    4*EPSILON*abs(x) + tol
  22.  *
  23.  * Algorithm
  24.  *    G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
  25.  *    computations. M., Mir, 1980, p.180 of the Russian edition
  26.  *
  27.  * The function makes use of the bissection procedure combined with
  28.  * the linear or quadratic inverse interpolation.
  29.  * At every step program operates three abscissae - a, b, and c.
  30.  *    b - the last and the best approximation to the root
  31.  *    a - the last but one approximation
  32.  *    c - the last but one or even earlier approximation such that
  33.  *        1) |f(b)| <= |f(c)|
  34.  *        2) f(b) and f(c) have opposite signs, i.e. b and c encompass
  35.  *           the root
  36.  * At every step Zeroin computes two new approximations, one by the 
  37.  * bissection procedure and the other one from interpolation (if a,b, and c
  38.  * are all different the quadratic interpolation is used, linear otherwise).
  39.  * If the latter (i.e. obtained by the interpolation) point looks
  40.  * reasonable (i.e. falls within the current interval [b,c] not close
  41.  * to the end points of the interval), the point is accepted as a new
  42.  * approximation to the root. Otherwise, the bissection result is used.
  43.  * Therefore, the range of uncertainty is guaranteed to be reduced at 
  44.  * least by the factor of 1.6
  45.  *
  46.  ************************************************************************
  47.  */
  48.  
  49. #include "math_num.h"
  50.  
  51. #pragma implementation "math_num.h"
  52.  
  53.  
  54. double zeroin(                // An estimate to the root
  55.     const double ax,        // Specify the interval the root
  56.     const double bx,        // to be sought in
  57.     double (*f)(const double x),    // Function under investigation
  58.     const double tol = EPSILON )    // Acceptable tolerance
  59. {
  60.   double a = ax, b = bx, c;        // Abscissae, see above
  61.   double fa;                // f(a)
  62.   double fb;                // f(b)
  63.   double fc;                // f(c)
  64.  
  65.   assure( tol > 0, "Tolerance must be positive");
  66.   assure( b > a, 
  67.      "Left end point of the interval should be strictly less than the "
  68.      "right one" );
  69.  
  70.   fa = (*f)(a);  fb = (*f)(b);
  71.   c = a;   fc = fa;
  72.  
  73.   for(;;)        // Main iteration loop
  74.   {
  75.     double prev_step = b-a;        // Step from the previous iteration
  76.    
  77.     if( fabs(fc) < fabs(fb) )
  78.     {                                 // Swap data so that b would be the
  79.       a = b;  b = c;  c = a;              // best approximation found so far
  80.       fa=fb;  fb=fc;  fc=fa;
  81.     }
  82.                     // Estimate the effective tolerance
  83.     const double tol_act = 2*EPSILON*fabs(b) + tol/2;
  84.     double new_step = (c-b)/2;        // Bissection step for this iteration
  85.  
  86.     if( fabs(new_step) <= tol_act || fb == 0 )
  87.       return b;                // Acceptable approximation is found
  88.  
  89.                 // Figuring out if the interpolation can be tried
  90.     if( fabs(prev_step) >= tol_act    // If prev_step was large enough
  91.     && fabs(fa) > fabs(fb) )    // and was in true direction,
  92.     {                    // Interpolatiom may be tried
  93.  
  94.       double p;                  // Interpolation step is calcu-
  95.       double q;                  // lated in the form p/q; divi-
  96.                       // sion operations is delayed
  97.                      // until the last moment
  98.       const double cb = c-b;
  99.  
  100.       if( a==c )            // If we've got only two distinct
  101.       {                    // points linear interpolation
  102.     register double t1 = fb/fa;    // can only be applied
  103.     p = cb*t1;
  104.     q = 1.0 - t1;
  105.       }
  106.       else                // Quadratic inverse interpolation
  107.       {
  108.     register double t1, t2;
  109.     q = fa/fc;  t1 = fb/fc;  t2 = fb/fa;
  110.     p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
  111.     q = (q-1.0) * (t1-1.0) * (t2-1.0);
  112.       }
  113.  
  114.       if( p > 0 )            // Formulas above computed new_step
  115.     q = -q;                // = p/q with wrong sign (on purpose).
  116.       else                // Correct this, but in such a way so
  117.     p = -p;                // that p would be positive
  118.       
  119.       if( p < (0.75*cb*q-fabs(tol_act*q)/2)    // If b+p/q falls in [b,c]
  120.      && p < fabs(prev_step*q/2) )    // and isn't too large
  121.     new_step = p/q;            // it is accepted
  122.                     // If p/q is too large then the
  123.                     // bissection procedure can
  124.                     // reduce [b,c] to a larger
  125.                     // extent
  126.     }
  127.  
  128.     if( fabs(new_step) < tol_act )    // Adjust the step to be not less
  129.       if( new_step > 0 )        // than the tolerance
  130.     new_step = tol_act;
  131.       else
  132.     new_step = -tol_act;
  133.  
  134.     a = b;  fa = fb;            // Save the previous approximation
  135.     b += new_step;  fb = (*f)(b);    // Do step to a new approximation
  136.     if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) )
  137.     {                             // Adjust c for it to have the sign
  138.       c = a;  fc = fa;                  // opposite to that of b
  139.     }
  140.   }
  141.  
  142. }
  143.